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We show that under certain circumstances the differences between the ab- 
sorption mean and Planck mean opacities can lead to multiple solutions for an 
CN \ LTE atmospheric structure. Since the absorption and Planck mean opacities 

are not expected to differ significantly in the usual case of radiative equilibrium, 
non-irradiated atmospheres, the most interesting situations where the effect may 
■^j- | play a role are strongly irradiated stars and planets, and also possibly structures 

where there is a significant deposition of mechanical energy, such as stellar chro- 
mospheres and accretion disks. We have presented an illustrative example of a 
strongly irradiated giant planet where the bifurcation effect is predicted to occur 
for a certain range of distances from the star. 
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1. Introduction 

There are many situations where the atmosphere of an object, a star or a planet, is 
irradiated by a companion star in such a way that this irradiation significantly influences 
its atmospheric structure. In the case of classical close binary stars, the effect exists, but is 
rarely dramatic because the effective temperatures differ by a factor of a few. On the other 
hand, the effect may be quite dramatic in the case of sub-stellar mass objects such as giant 
planets and brown dwarfs, irradiated by a solar-type star, in which case the ratio of their 
effective temperatures may reach a factor of 100 or more. We stress that we use the term 
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effective temperature as used in the theory of stellar atmospheres, namely as a measure of 
the total energy flux coming from the interior of the object. 

In the case of extrasolar giant planets (EGP), we recently have witnessed a significant 
increase in interest in predicting EGP spectra. This was motivated in part by two detected 
planetary transits, which are in principle able to provide direct information about the planet's 
atmosphere. Another motivation is the need to predict spectra of extrasolar planets to guide 
the design of future missions which aim at recording EGP spectra . 

We have recently computed a large set of model atmospheres and spectra of EGPs (Su- 
darsky, Burrows, & Hubeny 2003 - hereafter referred to as SBH). We have used a modification 
of our universal stellar atmosphere program TLUSTY (Hubeny 1988; Hubeny & Lanz 1995), 
called COOLTLUSTY. This variant does not compute opacities on the fly; instead it uses a 
pretabulated opacity as a function of wavelength, temperature, and density. The tables are 
computed using an updated version of the chemical equilibrium code of Burrows & Sharp 
(1999), which includes a prescription to account for the rainout of species in a gravitational 
field. 

Motivated by the recent discovery of the second transiting planet OGLE-TR-56 (Konacki 
et al. 2003; Sasselov 2003), which is believed to have a separation of a mere 0.0225 AU 
from its parent star, we tried to extend the SBH models to higher irradiations. However, 
COOLTLUSTY faced significant convergence problems. After various attempts and using 
various strategies, we have discovered that the convergence behavior is not the result of a bug 
in the program, or of an insufficiency in our numerical scheme, but is in fact a consequence of 
an existence of multiple solutions that may lead to disastrous effects for convergence. When 
studying the effect, we found that its roots are quite general, and are in fact applicable to 
other cases in the theory of stellar atmospheres. Therefore, we devote the present paper to 
explain the effect in general terms. 

2. Basic Equations 

The atmospheric structure is obtained by solving simultaneously the radiative transfer 
equation, the radiative equilibrium equation, and the hydrostatic equilibrium equation. For 
simplicity, we assume LTE. Since our basic aim here is to study atmospheres of irradiated 
giant planets, this approximation is a reasonable one, although it should be relaxed in the fu- 
ture, as was the LTE assumption eventually relaxed in modern studies of stellar atmospheres. 
Here, we present a brief overview of the basic equations. 
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The radiative transfer equation is written as 

where I Ufl is the specific intensity of radiation as a function of frequency, u, angle - described 
through the cosine of the angle of propagation with respect to the normal to the surface, 
/i, and the geometrical coordinate, taken here as the column mass m. The later is defined 
as dm = —pdz, where z is a geometrical distance (measured outwards), and p the mass 
density. The monochromatic optical depth is defined as dr u = Xu dm. Finally, S v is the 
source function, given in LTE by 

S v = ^-B v + ^-J v . (2) 

Xu Xv 

Here, k u is the true absorption coefficient, a u the scattering coefficient, and Xu = k u + cr u is 
the total extinction coefficient. All coefficients are per unit mass. 

The boundary conditions are provided through the diffusion approximation at the deep- 
est point, 

dB 

= ^(T(r max )) + fi , n > , (3) 

ar u 

and the upper boundary condition is 

/,„(0) = U, a*<0, (4) 

where 1^ is the specific intensity of the external irradiation. For simplicity, we assume 
isotropic irradiation, 1^ = J° xt . It is convenient to express the frequency-integrated irradi- 
ation intensity as 

poo 

jcxt = / jcxt dv = WB ( T ^ ? ( 5 ) 

Jo 

where is an effective temperature of the irradiating star (in case the irradiation source is 
not a star, is merely the characteristic temperature of the incoming radiation), and W a 
dilution factor. 

The moments of the specific intensity are defined as 

[J V ,H V ,K V ] = (1/2) j h,[l,p,p 2 ]dp. (6) 

The first moment of the transfer equation is written 

-7-^ = Xu {Ju - S v ) = K u (J v - B v ) , (7) 
dm 
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where the second equality follows from equation (2). Integrating over frequency we obtain 

^L = KjJ-K B B, (8) 

where Kj and k b are the absorption and Planck mean opacities, respectively, defined by 

_ j™K u J u dv 

"'J poo T , ) \y ) 

Jo J » dv 



and 



KB " / °° B v du ■ (10) 
These two opacities are usually assumed to be equal. This is an excellent approximation 
in the case of non-irradiated atmospheres, because J v ~ B v in optically thick regions, and 
J v [t < 1) oc (1/2)B v (t = 1) (the Eddington-Barbier relation), so J is proportional to B and 
the averaging over J and B leads to very similar results. However, here we maintain the 
distinction because the difference between kj and k b turns out to be crucial in the case of 
strongly irradiated atmospheres. 

The second moment of the transfer equation is 

dK v 



dm 

and integrating over frequency we obtain 



X,H„, (11) 



where 

_ j™XuH u dv 
Jo ^ rfz/ 

which is called the flux mean opacity. Notice that unlike the two previous opacities, which 
were averages of the true absorption coefficient (without the scattering term), the flux-mean 
opacity contains the total absorption coefficient. 

Finally, the radiative equilibrium equation is written 

POO 

/ Kv (J v -B v )du = 0, (14) 
Jo 

which can be rewritten, using the above defined mean opacities, as 

kjJ - k b B = 0. (15) 
Substituting (15) into (8), we obtain another form of the radiative equilibrium equation: 

*¥L = , or H = const = (a/Ait) T p 4 ff , (16) 
dm 

where a is the Stefan-Boltzmann constant. 
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3. Temperature Structure 
3.1. General 

The above equations are exact, given LTE, but are of course only formal because the 
mean opacities kj and xh are not known in advance; only k b is a known function of temper- 
ature and density. Nevertheless, assuming that they are known and equal to the Rosseland 
mean opacity, we may write a solution for temperature, following the classical textbook 
procedure, known by the name of "LTE-grey model atmospheres" (e.g. Mihalas 1978). A 
generalization of a classical model for the case of external irradiation is given by Hummer 
(1982), and for the case of accretion disks and unequal mean opacities by Hubeny (1990). 

The procedure is as follows. From equation (15) we have B = (k,j/k,b)J, which allows 
us to express T through J using the well-known relation B = (a/njT 4 . To determine J, we 
use the solution for the second moment of the transfer equation K(th) = Hth + K(0) = 
(a I in) T e g th + K(0), where th is the optical depth associated with the flux-mean opacity, 
and express the moment K in terms J using the Eddington factor, fx = Kj J. Similarly, 
we express the surface flux through the second Eddington factor, f H = iJ(0)/J(0), and we 
end up with (see also Hubeny 1990) 



T 4 Cff ^ 



1 

-th + 



+ ^WT*. (17) 
Kb 



3f K 3/ H . 

Again, this solution is exact within LTE. The usual LTE-gray model consists in assuming 
that all the mean opacities are equal to the Rosseland mean opacity. Moreover, if one adopts 
the Eddington approximation (f K — 1/3; fn — l/VS), then one obtains a simple expression 

^ 3 



rpA 
l 1 ^ 



(r+l/Vt)+WTt. (18) 



We will consider the most interesting case, namely that of strong irradiation, defined by 
WT 4 3> T e 4 ff . In this case, the second term in brackets is negligible, and we may define 
a penetration depth as the optical depth where the usual thermal part (oc T c 4 ff ) and the 
irradiation part (oc WT*) are nearly equal, viz. 

r pen » W (J^J . (19) 

The behavior of the local temperature in the case of the strictly gray model is very simple - 
it is essentially constant, T = T = W 1 ^^ for r < r pcn , and follows the usual distribution 
T oc r^Tefj in deep layers, r > r pcn . We stress that while the strictly gray model exhibits 
an essentially isothermal structure down to r ^ r pcn , such will not be the case for a nongray 
model, as we shall show in detail in the next sections. 
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3.2. Surface Layers 

In the general case, we have to retain the ratio of the absorption and Planck mean 
opacities. In the irradiation-dominated layers (t# < r pcn ), the temperature is given by 

r = 7 r /4 r i , (20) 

where 

^^{kj/kb) 1 ' 4 . (21) 

As stated before, 7 « 1 in the case of no or weakly irradiated atmospheres. However, in 
the case of strong irradiation, 7 may differ significantly from unity, and, moreover, may 
be a strong function of temperature, and to a lesser extent density. This is easily seen by 
noticing that in optically thin regions, the local mean intensity is essentially equal to twice 
the irradiation intensity, because the incoming intensity is equal to the irradiation intensity, 
and the outgoing intensity is roughly equal to it as well. The reason is that in order to 
conserve the total flux when it is much smaller than the partial flux in the inward or the 
outward direction, both fluxes should be almost equal, as are the individual specific and 
mean intensities. More specifically, 



H = H out -H cxt = / I^fidfidu- / / I u ^fidfxdu. (22) 
Jo Jo Jo Jo 

If H <C H cxt , then we must have I vll dv J °° I^^du for all angles, and thus J v « 

2jext ^ 2 W x ^B v iJ^^). We may then write the absorption mean opacity as a function of T 

and T* as 

(TT , fKy(T)WB v (T.)du _ J K v (T)B v (T*)dv 
M ' * j ~ fWB v (T.)dv ~ fB u (T*)dv ■ [I6) 
The dilution factor cancels out, and only the spectral distribution of the irradiation intensity 
matters. 

In other words, the absorption mean is an average of the opacity weighted by the Planck 
function corresponding to the effective temperature of the source of irradiation, while the 
Planck mean is an average of the opacity weighted by the Planck function corresponding to 
the local temperature. Obviously, they can be quite different. If the monochromatic opacity 
differs significantly in the region where B V {T) and B U (T*) have their local maxima, the 
resulting kj and Kb will differ substantially. If, moreover, and this is the crucial point, the 
monochromatic opacity depends sensitively on temperature, the opacity ratio 7 may have a 
complex, and generally non-monotonic, dependence on temperature. 

The local temperature in the upper layers is given, in view of equation (20), by the 
expression 

T/T = 7 (T) . (24) 
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It is now clear that if 7 exhibits a strongly non-monotonic behavior in the vicinity of T , for 
instance if it has a pronounced minimum or maximum there, equation (24) may have two or 
even more solutions] 



3.3. Deep Layers 

The bifurcation behavior is not limited to the surface layers, but may continue to large 
optical depths. This might seem surprising at first sight because from Equation (17) one 
may expect that once r > 1, the mean opacities kj and Kb become roughly equal, and 
the Eddington factor is fx ~ 1/3, so that T = T all the way till r pen . We recall that the 
penetration depth may be quite large; for instance for the case T e g = 75 K, T* = 6000 K, 
and W = 2.2 x 10~ 3 (the case studied in detail in the next Section), r pcn = 9 x 10 4 . 

However, the behavior of the local temperature may be more complicated. This is 
linked to another interesting inequality of mean opacities that is usually taken for granted, 
namely the flux mean opacity xh and the Rosseland mean opacity Xross- The Rosseland 
mean opacity is in fact defined in such a way that it is equal to the flux mean opacity in the 
diffusion approximation. Indeed, in this approximation, 

\dB v 11 dB u 1 1 dB u olT 

J J r^j 

3 dr v 3 Xv dm 3 Xu dT dm ' 

and therefore 

Jo°° IFF*" 
v ~ J{ J til _ 

where the latter equality represents the definition of the Rosseland mean opacity. 



In the case of strong irradiation, an interesting, and fundamentally different, situation 
appears. Since the net flux is very small, the total flux in the fi > (outgoing) hemisphere 
is roughly equal to that in the < (incoming) hemisphere (see above). Because the 
monochromatic opacity varies strongly with frequency, there are frequency regions where 
the net monochromatic flux is positive, and regions where it is negative. The flux-mean 
opacity close to the surface may attain large values, either positive or negative, depending 
upon whether Xv weighs the positive or negative net flux regions more. 

Going to deeper layers, the net monochromatic flux decreases because the radiation 
field becomes more isotropic for all frequencies, so that the flux mean opacity decreases. 
Consequently, the corresponding flux-mean optical depth increments Ath ~ xu^m become 
very small, and th will exhibit a plateau where it remains essentially constant with m. 
Finally, in the regions where all the influence of the external irradiation dies out, the usual 
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diffusion approximation sets in, and the flux-mean optical depth becomes essentially equal 
to the Rosseland optical depth. 

This means that for strong irradiation, the temperature in the deep layers should exhibit 
a plateau with 

3 

Tplatcau ~ | T* s f H , (25) 

where f# is the plateau value of the flux-mean optical depth. The value of f H and thus 
7pi a teau depends upon the exact form of the monochromatic opacity; a rough estimate is 
provided by setting f H pz r pen , in which case T platC a U ~ (7/4) 1 / 4 W 1 / 4 T„, 1.15 T . 



4. Example of an Extrasolar Giant Planet Irradiated by a Solar-type Star 

4.1. Opacities 

Using partition functions, LTE level densities, stimulated emission corrections, and 
broadening algorithms, we generated opacity tables from numerous available line lists. For 
gaseous U2O, Partridge and Schwenke (1997) have calculated the strengths of more than 
3 x 10 8 lines. We used a subset of their 4 x 10 7 strongest lines. For various other molecular 
species (e.g., NH 3 , PH 3 , H 2 S and CO), we used the HITRAN (Rothman et al. 1992,1998) 
and GEISA (Husson et al. 1994) databases, augmented with additional lines from theoret- 
ical calculations and measurements (Tyuterev et al. 1994; Goorvitch 1994; Tipping 1990; 
Wattson and Rothman 1992). For methane, shortward of ~1.0 /iinwe used the Karkoschka 
(1994) opacities and between ~1.0 /inland 1.58 /imwe used the Strong et al. (1993) opac- 
ities. Longward of 1.58 /im, we used the Dijon methane database (Borysow et al. 2003) 
that includes the hot bands. Opacity due to Collision-Induced Absorption (CIA) by H 2 
and helium is taken from Borysow and Frommhold (1990), Zheng and Borysow (1995), and 
Borysow, J0rgensen, and Zheng (1997), as updated in Borysow (2002). 

FeH and CrH opacities were taken from Dulick et al. (2003) and Burrows et al. (2002), 
respectively. The TiO line lists for its nine major electronic systems were taken from the 
Schwenke ab initio calculations, as modified by Allard, Hauschildt, and Schwenke (2000). 
These data include lines due to isotopically substituted molecules 46 Ti 16 0, 47 Ti 16 0, 49 Ti 16 
and 50 Ti 16 O relative to the most abundant isotopic form 48 Ti 16 0. For VO, we used the line 
list provided by Plez (1999). Because 51 V is by far vanadium's most abundant isotope, the 
lines of isotopically substituted molecules are not necessary. The line lists and strengths for 
the neutral alkali elements (Li, Na, K, Rb and Cs) were obtained from the Vienna Atomic 
Line Data Base (Piskunov et al. 1995). The general line shape theory of Burrows, Marley, 
and Sharp (2000) was used for Na and K, while those of Nefedov et al. (1999) and Dimitrijevic 
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and Peach (1990) were used for the other alkali metals. 

In practice, we pre-calculate a large opacity table in T/p/frequency space in which 
we interpolate during the iteration of the atmosphere/spectral model to convergence. The 
table contains 30,000 frequency points from 0.3 to 300 pm, uniformly spaced in steps of 
1 cm -1 . Generally, 5000 geometrically-spaced frequency points are used in the radiative 
transfer model. In other words, our approach belongs to the category of "opacity sampling" 
methods for treating line blanketing. 

We consider two different opacity tables; one which includes opacity due to TiO and VO, 
and the other one without these molecules. In the case without TiO/VO, these species were 
removed both from the opacity table as well as from the state equation, so these species are 
considered as completely absent everywhere in the atmosphere. These two different opacity 
tables suggest themselves because for an irradiated atmosphere we obtain high temperatures 
in the low pressure outer atmosphere. Importantly, temperature inversions are very possible. 
However, thermochemical equilibrium calculations with rainout (Burrows & Sharp 1999; 
Lodders & Fegley 2002) suggest that TiO and VO would exist at such high-T/low-P points. 
Since the presence of TiO and VO is quite natural at high-T/high-P points, one confronts a 
situation in which there are two regions rich in TiO/VO, separated by a middle region that 
is TiO/VO free. In a gravitational field with a monotonic pressure profile, this gap could 
act like a cold trap in which the TiO/VO that is transported by molecular or eddy diffusion 
from the upper low-P region into the intermediate cooler region, would condense and settle 
out, thereby depleting the upper low-P TiO/VO-rich region. This would eventually leave no 
TiO/VO at altitude to provide the significant absorption that could lead to a bifurcation. 
However, the cold-trap effect can not be anticipated in abundance tables that are functions 
of temperature and pressure alone. These tables are not cognizant of gravity, nor are they 
aware of the global T/P profile. Whether this cold-trap effect in fact happens is not yet 
clear. Furthermore, the same cold-trap effect might obtain for other EGP species (e.g. Fe). 
Hence, we use the presence or absence of TiO/VO at high-T/low-P points and the associated 
ambiguity as a means to explore the real mathematical bifurcation effect we have identified 
in this paper and leave to a future work the study of the true viability of cold traps in 
irradiated EGP atmospheres. 



4.2. Results 

We consider an intrinsically cold giant planet, with T c g = 75 K, logg = 3, irradiated by 
a GOV star with a spectral energy distribution taken for simplicity to be the corresponding 
Kurucz (1994) model atmosphere. We consider the class V, "roaster" situation, as defined 
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by SBH. 

The opacity ratio index 7 is plotted in Figure 1 as a function of temperature, for three 
densities characteristic to the outer layers of an atmosphere - p = 10~ 12 , 10 -11 , and 10 -10 g 
cm -3 . The upper panel shows the results for the opacities including TiO/VO, while the lower 
panel shows those without TiO/VO. It is clearly seen that in the former case the ratio k,j/k b , 
and thus 7, exhibits a sharp peak around 1500 K; the dependence on density is not strong. 
The ratios T/T for T = 450 K (weak irradiation), 1250 K (medium strong irradiation), 
and 2600 K (very strong irradiation) are overplotted with dotted lines. The latter roughly 
corresponds to the case of the recently discovered transiting planet OGLE-TR-56. According 
to equation (24), the temperature at the surface layers is found at the intersection of the 
curves of j(T) and T/T Q ; thus, we see that for low and for high irradiation there are unique 
solutions (corresponding toTsi 250 K and T m 2600 — 30007^, but that in the intermediate 
cases there are indeed multiple solutions. 

The opacity table without TiO/VO does not exhibit a peak around 1500 K, and as it is 
clearly seen from the lower panel of Fig. 1, one obtains unique solutions for the temperature 
for all reasonable irradiations. Mathematically, the case of truly extreme irradiation (e.g. 
with T around 4700 K) might again lead to multiple solutions 4 because the line T/T now 
almost coincides with 7(T). However, we do not consider such a case here because then a 
host of new phenomena which we do not address in our present model, such as departures 
from LTE, dynamical effects, etc., would likely become crucial. 

To explain the behavior of the kj over Kb ratio, we plot in Figure 2 two examples of 
the monochromatic opacity, for T = 1600 K (close to the maximum of 7), and T = 520 
K (the region near the minimum of 7). In both cases p = 10~ n g cm~ 3 . We also plot 
normalized Planck functions for the two nominal temperatures, 1600 and 520 K (dotted 
lines), and for T = T* = 6000 K (dashed line). The latter weighs most heavily the optical 
region around v pa 4 — 5 x 10 14 Hz, while the local Planck functions put the weight at 
much lower frequencies (around 10 14 and 3 x 10 13 Hz, respectively). In the case of the higher 
temperature, T = 1600 K, the monochromatic opacity exhibits a maximum just in the region 
where the Planck function corresponding to T* is largest (because of TiO and VO); therefore, 
Kj > k,b- For the lower temperature, T = 520 K, the opacity in the optical range is lower by 
orders of magnitude, so kj is now significantly lower. Since the Planck function for T = 520 
K emphasizes frequencies below 7 x 10 13 Hz, where the monochromatic opacity is very large, 
we obtain kj < Kb- 

There are two other features seen in Figure 2 that are worth noticing. First, there is 
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a large increase of 7 for temperatures around and below 100 K. This is a consequence of 
low opacity at low frequencies. However, the opacity sources are quite uncertain in this 
regime. For studies of extremely cool objects the opacity table should be checked and 
possibly upgraded, but studying these objects is not our objective here. Second, in the case 
of no TiO/VO opacity (lower panel), it appears that 7 nearly coincides with T/T*. Taken 
at the face value, this would imply that the ratio J K v B v {T^)dv / J n v B v (T)dv is very close 
to unity. 5 However, we do not see any compelling reason why this should be generally valid, 
so we conclude that such a behavior is just a coincidence. 

In Figure 3 we demonstrate the behavior of the flux-mean opacity discussed in §3.3. 
The high-surface-T branch exhibits a negative flux mean optical depth at the surface layers 
because the opacity is high in the optical region where the next flux is negative. Conse- 
quently, the flux-mean opacity is negative there, and thus the flux-mean opacity (and the 
corresponding optical depth) on the plateau is lower than in the case of the low-surface-T 
branch. 

We have also computed exact LTE models as described in SBH. We have used the stellar 
atmosphere code COOLTLUSTY, which is a variant of the universal program TLUSTY 
(Hubeny 1988; Hubeny & Lanz 1995). We stress that the emergent flux is not computed 
by a separate program; instead it is computed already by COOLTLUSTY. This also means 
that we do not use different opacity tables for computing the atmospheric structure and 
for computing the spectra. We have also tested the sensitivity of the computed model on 
different opacity samplings; we have resampled our original frequency grid of 5000 points to 
lower numbers of frequencies, and even for 300 points we found very little differences in the 
temperature structure (of the order of a few K). In some cases it was found advantageous 
to converge first a model with 300 frequencies, and using this as an input to reconverge a 
model with the full grid of 5000 frequencies in the next step. 

The temperature structure is displayed in Figure 4, which demonstrates that the simpli- 
fied description put forward above faithfully reflects the true behavior of temperature. Full 
lines represent models computed for the opacity table without TiO/VO, while the dashed 
and dot-dashed lines represent models with TiO/VO. The models without TiO/VO always 
converged to a unique solution, moreover with temperature monotonically decreasing out- 
ward, as we may have expected from the above discussion. 

The behavior of the models with TiO/VO is more interesting. The high-irradiation 
model, the top dashed line, indeed converged to a unique solution. The surface temperature 
is 2600 K, in excellent agreement with the value expected from Figure 1 (the density at the 
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uppermost point is about 3 x 1CT 12 g cm -3 ). Even the initial increase of temperature in 
the upper layers when going inward is explained by Figure 1, since for progressively higher 
densities the solution of T/T = 7(T) occurs at higher temperatures. The low-surface-T 
solution is very close to the no-TiO/VO model, which once again demonstrates that the 
bifurcation is caused by the TiO/VO opacity. We have also examined depth-dependent 
concentrations of TiO and VO. Indeed, the mole fraction of TiO drops precipitously from 
3 x 1CT 7 to essentially zero at m ^ 1.6 x 1CT 2 g cm~ 3 (that is at pressure « 1.6 x 1CT 5 
bars and temperature T m 1500 K) for the high-surface-T model; while is is not present in 
the low-surface-T model. Behavior of the VO concentration is analogous. 

The low-irradiation model also converged to a unique solution; moreover, its surface 
temperature of about 250 K is in excellent agreement with the value one may expect from 
Figure 1. Again, both models, without and with TiO/VO agree very well. 

The two middle curves, the dashed and dot-dashed curves, correspond to the same irra- 
diation, and, thus, represent the bifurcation predicted by our analysis. The model with high 
surface temperature was converged starting from the extremely irradiated model (through 
several intermediate steps), while the model with the low surface temperature was converged 
starting from scratch, i.e. from a usual LTE-gray model in which all mean opacities are cre- 
ated equal. The high surface temperature at about 2200 K essentially coincides with the 
high-temperature intersection of about 2200 K on Figure 1. The low-surface-T solution 
exhibits a somewhat higher temperature than would follow from the low-T intersection in 
Figure 1 - 650 versus 550 K, but in view of all the approximations involved such an agree- 
ment is still very good. We did not numerically find the third possible solution with a surface 
T of about 1300 - 1400 K. 

We stress that the numerical solutions were not contrived in any way to lead to the 
predicted bifurcation. In fact, we first discovered the effect purely numerically - we obtained 
two different, yet perfectly well converged solutions, depending on the input model. We have 
continued the convergence of both solution until the maximum relative change of temperature 
and density decreased below 10~ 5 in all depth points. 

Several features of the overall accuracy of the model are worth mentioning. In fact, 
the low value of maximum relative change of temperature and density is not a satisfactory 
criterion of the convergence. One should also check the accuracy of the computed flux 
gradient and the value of the net flux. That is, for the former case we have to examine the 
quantity e grad = (kjJ—k b B)/kbB 1 and for the latter case the quantity e flux = (H—H )/H ), 
where H = (a/4ir)T^ s is the net flux. For the models presented in Figure 4, e gra d is very 
small, typically 10~ 10 or less. The conservation of the total net flux is much harder to achieve; 
the bifurcated models have typically efl ux ~ few x 10~ 4 near the surface, while efl ux reaches a 
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few xl0~ 2 for column mass around m « 10 2 g cm -2 (that is, a few per cent), and drops to 
less than 10~ 5 at m ^ 10 4 g cm -2 and deeper. The error in the net total flux of several per 
cent might seem large, but one has to bear in mind that the net flux is a subtraction of two 
large quantities, H ont — H m , and the actual numerical error (H out — H U1 — H )/ H out is of the 
order of 10~ 5 or less. A dramatic demonstration of this fact is that a model, which is almost 
converged (maximum relative change in T and p was less than 10 -3 ), exhibits a seemingly 
ridiculous value of e flux « 10 5 at upper layers m < 10 2 g cur 2 (i.e. an error in the total 
net flux of some 10 7 % !), yet the actual temperature difference between this and the fully 
converged model is it most about 2-3 K, and the predicted flux in both cases is completely 
indistinguishable on a plot! This again shows that in the case of strong irradiation it is the 
value of (H out — H m — H )/H out (or e gra d) which is critical for practical convergence, not a 
much more stringent criterion of the error in the total net flux, efl ux . 

The difference in the temperature structure is of course reflected in the emergent spec- 
trum. In Figure 5 we show the emergent flux corresponding to the two solutions of the 
intermediate-irradiation (distance 0.08 AU) model, computed with TiO/VO. The models 
differ dramatically in the optical and in near-IR region because of the dramatically different 
surface temperature and the corresponding increase of the TiO/VO opacity. Since the low- 
surface- T branch exhibits a higher temperature at the plateau (see Fig. 4), the flux at the 
low-opacity regions is higher. 

Finally, Fig. 6 displays the emergent flux for the two models for the highest irradiation, 
namely for the distance 0.0225 AU, which corresponds to OGLE-TR-56. The thick line is a 
model without TiO/VO, and the thin line with TiO/VO. There is no bifurcation here since 
both opacity tables led to unique solutions, however quite different ones. The corresponding 
spectra are thus significantly different as well. 

5. Discussion and Conclusions 

We have demonstrated that under certain circumstances the differences between the 
absorption mean and the Planck mean opacities can lead to multiple solutions for the at- 
mospheric structure. This result is quite robust, and in fact represents a so far overlooked 
general result of elementary stellar atmospheres theory. Since we do not expect that the 
absorption and Planck mean opacities will differ significantly in the usual case of radiative 
equilibrium of non-irradiated atmospheres, the most interesting situations where the effect 
may play role are strongly irradiated stars and planets, and also, possibly, structures where 
there is a significant deposition of mechanical energy, such as stellar chromospheres and 
accretion disks. We are not concerned with these objects here, but would like to stimulate 
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further study of possible bifurcations for these situations. We note that Nolan & Lunine 
(1988) also found a bimodal behavior of the atmosphere of Triton which is linked to external 
irradiation. 

From the physical point of view we may interpret the absorption mean opacity as the 
global absorption efficiency, and the Planck mean opacity as the global emission efficiency, 
of the medium. The integrated mean intensity J acts as a total absorption pool, while 
the integrated Planck function B as a total thermal pool. The radiative equilibrium, that 
stipulates that the total radiation energy absorbed on the spot is equal to the total energy 
emitted on the same spot, therefore acts as a thermostat: The radiative equilibrium sets the 
local temperature in such a way that KjJ = kbB. Here J is determined by the radiative 
transfer, and the local temperature follows as B = (cr/7r)T 4 = (kj/k b )J. In the case of 
strong irradiation, the weighting that determines Kj is dominated by the incoming radiation. 
If the monochromatic opacity in the region around the peak of the external irradiation is 
very sensitive to temperature, for instance if it is very high for high temperatures and very 
low for low temperatures, the "radiative-equilibrium thermostat" may find two solutions - 
either a high T together with high kj, or a low T with low kj] in both cases the radiative 
equilibrium KjJ = k b B is satisfied. 

However, the applicability of this effect to real objects depends critically on a degree of 
reliability of adopted opacities. For high-temperature conditions where the gaseous opacities 
dominate, the situation is relatively well under control, thanks to recent enormous progress 
in computing atomic data (Opacity Project, Iron Project, OPAL, and others; for a recent 
review see, e.g., Nahar 2003). 

In the case of irradiated giant planets or other substellar-mass objects like L and T 
dwarfs, the opacity for low density (p ~ 10~ 12 — 10~ 9 g cm" 3 ) and high temperature (T > 2000 
K) is still relatively uncertain. At the present stage, we cannot be sure that the bifurcation 
effect really occurs, but we have shown that this is certainly a very real possibility which 
should be taken into account when constructing atmospheric models. 

The authors are pleased to thank Bill Hubbard, Thierry Lanz, Jonathan Lunine, Jim 
Liebert, Christopher Sharp, and Drew Milsom for fruitful conversations and help during the 
course of this work, as well as NASA for its financial support via grants NAG5-10760 and 
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Fig. 1. — Opacity ratio index 7 = (kj/k b ) 1 I a plotted as a function of temperature. The 
upper panel shows the results for the opacities including TiO/VO, while the lower panel 
shows those without TiO/VO. The absorption mean opacity was evaluated using equation 
(23) with the stellar effective temperature T* = 6000 K. The full line corresponds to a 
density p of 10~ 12 g cm" 3 ; the dashed line to p = 10~ n g cm~ 3 ; and the dot-dashed line 
to p = 10~ 10 g cm -3 . The dotted lines represent the quantity T/T , with T corresponding 
to three different distances from the star. From top to bottom, To = 450 K, (0.65 AU), 
T = 1250 K (0.08 AU), and T = 2600 K, (0.0225 AU). The intersections of the individual 
dotted lines with the j(T) curve represent possible solutions for the temperature in the 
upper layers of an irradiated atmosphere. In the case of the opacity table with TiO/VO, 
there are unique solutions for low and for high irradiation (corresponding toT« 250 K and 
T 2600 — 3000fT), while in the intermediate cases there are indeed multiple solutions. 
The opacity table without TiO/VO does not exhibit a peak around 1500 K, so one obtains 
unique solutions for the temperature for all irradiations. 
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Fig. 2. — Monochromatic opacity (in cm 2 g -1 ) as a function of frequency for density p = 
10~ n g cm' 3 and for two temperatures, T = 1600 K (solid line, with a maximum around 
v ~ 3 — 6 x 10 14 Hz), and T = 520 K (bold solid line); together with normalized Planck 
functions for T = T* (dashed line), and for the two nominal temperatures (1600 and 520 K 
- dotted lines). A large monochromatic opacity for T = 1600 K in the optical region where 
B V {T = T*) has a maximum explains why kj > k b for this temperature. 
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Fig. 3. — Flux-mean optical depth as a function of Rosseland mean optical depth for the 
two solutions of the intermediate-irradiation (distance 0.08 AU) model, computed with the 
TiO/VO opacity. The heavy line represents the low-surface-T branch; the thin line the high- 
surface-T branch. The solid line represents the positive values of t h , while the dashed line 
the negative values. 
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Fig. 4. — Temperature as a function of column mass. Solid lines represent models computed 
for the opacity table without TiO/VO, while the dashed and dot-dashed lines represent 
models with TiO/VO. The curves represent, from top to bottom, three distances from the 
central star. The top lines corresponds to 0.0225 AU, the middle lines to 0.08 AU, and 
the bottom lines to 0.65 AU. Numerical values of the dilution factor are W = 2.8 x 10~ 2 , 
2.2 x 10~ 3 , and 3.5 x 10~ 5 . The models for the lowest irradiation (W = 3.5 x 10~ 5 ) computed 
with and without TiO/VO are not distinguishable on the plot. The middle dashed and the 
dot-dashed lines represent the two solutions for the same distance, and thus illustrate the 
bifurcation discussed in the text. 
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Fig. 5. — Emergent flux of a function of wavelength for the two solutions of the intermediate- 
irradiation (distance 0.08 AU) model, computed with TiO/VO. Heavy line represents the 
low-surface-T branch; the thin line the high- surf ace-T branch. The model flux computed 
without TiO/VO is very close to the low-surface-T branch, and for clarity is not displayed. 
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Fig. 6. — Emergent flux of a function of wavelength for the high-irradiation (distance 0.0225 
AU) model, computed without TiO/VO (thick line) and with TiO/VO (thin line). 



